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Rayleigh  waves  normally  incident  upon  2-D  shallow  heterogeneity  are  simulated  by  the 
linear  finite-difference  method  to  study  attenuation,  transmission,  and  reflection  of  Rayleigh 
waves  and  to  measure  the  Rayleigh-to-P  and  -SV  body  wave  conversion.  Transmission, 
reflection,  and  scattering  depend  on  the  depth,  average  scale  size  of  the  heterogeneity  and  the 
amplitude  of  the  spatial  fluctuation  of  velocity.  As  expected,  larger  spatial  variation  in  velocity 
attenuates  Rayleigh  waves  more  than  smooth  media,  and  the  attenuation  is  roughly  proportional 
to  the  variance  of  the  velocity  fluctuation.  The  attenuation  and  scattering  due  to  shallow 
heterogeneity  is  weaker  than  attenuation  due  to  moderately  rough  topography. 

Scattered  body  wave  energy  is  studied  as  a  function  of  frequency,  scattering  angle,  and 
wave  type  (P  or  SV).  Attenuation  of  Rayleigh  waves  by  scattering  from  2-D  shallow  velocity 
heterogeneity  is  dominated  by  conversion  to  body  waves  and  in  particular  SV  energy.  Low 
frequency  P  and  SV  energy  is  scattered  in  a  backwards  direction,  and  high  frequency  P  and  SV 
energy  is  scattered  in  a  forward  direction. 

As  with  scattering  from  rough  topography,  much  of  the  converted  SV  energy  will  be 
trapped  in  the  crustal  waveguide  at  Lg  phase  velocities.  Therefore,  Rayleigh  (Rg)  to  SV  body 
wave  conversion  by  shallow  heterogeneity  and  topography  should  contribute  to  the  formation 
of  Lg  by  explosions,  quarry  blasts,  and  shallow  earthquakes. 

A  comparison  is  made  with  results  for  P-coda  from  Greenfield  (1971).  The  comparison 
indicates  that  self-similar  and  Gaussian  models  could  be  derived  with  rms  velocity  variations 
between  7  and  15%  in  the  upper  3  km  of  the  crust  that  would  produce  the  observed  P-coda/P 
power  levels  observed  by  Greenfield  (1971). 
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INTRODUCTION 

A  significant  fraction  of  the  elastic  energy  generated  by  a  shallow  explosion  is  contained 
in  short-period  Rayleigh  waves.  These  short-period  Rayleigh  waves  are  the  largest  phase  on 
the  near-source  seismogram.  Although  short-period  Rayleigh  or  Rg  are  commonly  seen  on 
regional  seismograms  from  shallow  sources  such  as  quarry  blasts,  they  are  attenuated  rather 
rapidly  with  spatial  Q’s  at  1  Hz  in  the  range  of  50  to  200  (Kafka,  1987;  and  personal  commun¬ 
ication). 

Attenuation  can  be  due  to  intrinsic  shear  attenuation  in  the  near  surface,  or  it  can  be 
accomplished  by  scattering  from  shallow  heterogeneity  and  topography.  If  the  energy  is  scat¬ 
tered  then  it  must  be  observed  as  coda  in  the  near-source,  regional,  and/or  teleseismic  seismo¬ 
gram.  Scattering  of  Rayleigh  waves  has  been  hypothesized  to  produce  teleseismic  P-wave 
coda  from  Novaya  Zemlya  explosions  by  Greenfield  (1971).  Gupta  et  at.  (1985)  and  Gupta  and 
Blandford  (1987)  propose  a  model  that  suggests  that  P-wave  coda  is  largely  composed  of  scat¬ 
tered  Rayleigh  waves  and  therefore  coda  should  contain  information  about  the  depth  of  burial 
of  an  explosion.  They  suggest  that  P-wave  coda  by  nature  of  its  averaging  properties  gives  a 
more  precise  measure  of  yield  than  the  direct  P  wave.  However,  since  we  have  no  accepted 
model  for  the  generation  of  P-coda  by  explosions,  such  a  relationship  must  be  established 
empirically  for  each  testing  area. 

In  support  of  the  Gupta  era/.  (1985)  model  for  Rayleigh  to  P-wave  coda  generation,  Lay 
(1987a)  observed  frequency  dependence  of  the  temporal  energy  centroid  from  Pahute  Mesa. 
Recently,  Lay  (1987b)  has  attempted  to  image  scattered  Rayleigh  waves  at  Pahute  Mesa 
assuming  that  short-period  (I  Hz)  teleseismic  P-wave  coda  is  produced  by  conversion  from 
Rayleigh  waves  initially  generated  by  the  explosions.  He  finds  that  1  Hz  coda  does  have  some 
spatial  patterns  reminiscent  of  the  focusing-defocusing  patterns  of  direct  P-waves  at  Pahute 
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Mesa.  Cormier  (1987)  puts  forth  a  model  for  coda  generation  that  hypothesizes  that  coda  only 
partially  averages  the  focusing-defocusing  near  the  source  and  near  the  receiver. 

McLaughlin  et  al.  (1983)  examined  coda  of  the  near-source  seismogram  at  Pahute  Mesa. 
Although  they  found  that  the  major  transverse  motion  observed  in  the  near-field  appeared  to 
arise  directly  from  the  explosion,  much  of  the  transverse  coda  energy  was  characteristic  of 
scattered  phases,  and  the  important  scale  lengths  of  shallow  heterogeneity  at  Pahute  Mesa  were 
in  the  0.5  to  1  km  range.  Such  scale  lengths  should  be  optimal  for  scattering  the  1  to  2  Hz 
Rayleigh  waves  at  Pahute  Mesa.  However,  since  the  mesa  has  both  topographic  relief  as  well 
as  considerable  geologic  complexity,  it  is  uncertain  at  this  time  as  to  what  structures  are 
responsible  for  the  teleseismic  P-coda. 

In  a  related  matter.  Key  (1967  and  1968)  observed  P  to  Rayleigh  scattering  by  specific 
isolated  topographic  features  in  the  vicinity  of  the  EKA  array.  Since  this  is  the  reciprocal 
problem  of  Rayleigh  to  P  wave  scattering,  one  can  consider  this  as  additional  motivation  for 
studying  the  phenomenon.  Also  using  array  techniques.  Dainty  (1986)  observed  that  much  of 
the  locally  generated  P-coda  at  NORSAR  is  characterized  by  energy  with  Rayleigh  wave  phase 
velocities.  Sometimes  the  NORSAR  coda  could  be  seen  to  arrive  from  a  specific  direction,  but 
it  was  often  isotropic.  Der  et  al.  (1984)  concluded  using  deconvolution  techniques  at  arrays 
that  the  vertical  P-wave  coda  consisted  of  approximately  50%  locally  generated  response  to 
plane  waves  arriving  near  the  receiver  and  of  approximately  50%  that  arrived  from  the  source 
region  as  plane  waves  across  the  array  at  the  same  slowness  as  the  direct  P  wave. 

The  fact  that  explosion  transverse  Lg  is  on  average  just  as  large  as  the  radial  or  vertical 
component  of  Lg  serves  to  highlight  our  lack  of  understanding  of  the  scattering  problem. 
Gupta  and  Blandford  (1983)  proposed  that  P-to-SV  and  SV-to-SH  scattering  was  responsible 
for  the  growth  of  the  transverse  Lg  with  distance,  and  Baumgardt  (1985)  has  proposed  that  P- 
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coda  at  NORSAR  can  be  observed  scattered  from  "Lg"  crossing  the  Ural  mountains.  Both  of 
these  problems  have  in  common  the  need  to  quantify  the  conversion  efficiency  of  modal  waves 
to  body  waves  and  vice  versa. 

The  generation  of  SV-Lg  by  explosions  is  an  important  problem  not  only  because  it  is 
poorly  understood  but  because  Lg  is  commonly  used  to  estimate  the  yield  of  explosions  (as  in 
Nuttli,  1986a,  1986b,  1987a,  and  1987b).  The  attenuation  correction  for  the  Lg  amplitude  is 
derived  by  an  analysis  of  the  Lg  coda  as  if  it  were  backscattered  Lg  modal  waves.  This  Lg 
coda  analysis  is  questioned  by  Der  et  al.  (1987).  They  claim  that  the  Lg  coda  is  converted  fun¬ 
damental  Rayleigh  and  not  proper  "Lg". 

Aki  and  Chouet  (1975)  proposed  models  for  the  generation  of  coda  as  the  backscattering 
of  S  waves  and/or  surface  waves.  Their  coda  model  of  S-wave  single  backscattering  has 
become  a  frequently  used  paradigm;  however,  their  observations  that  coda  levels  were  higher  at 
stations  located  on  alluvium  and  in  valleys  suggest  that  much  of  the  coda  is  surface- wave 
related.  Phillips  and  Aki  (1986)  found  additional  evidence  for  surface  waves  at  a  number  of 
such  stations.  Estimates  for  generation  of  surface  waves  by  incident  SV  waves  and  back- 
scattering  cross-sections  for  Rayleigh  waves  upon  topography  and  shallow  heterogeneity  would 
be  useful  in  the  modeling  of  coda.  Since  coda  Q  inferred  from  the  coda  of  seismograms  is 
often  interpreted  as  the  crustal  shear-wave  Q,  it  is  important  that  the  correct  model  for  the 
creation  and  propagation  of  coda  be  determined. 

Frankel  and  Wennenberg  (1987)  propose  a  model  for  coda  that  stresses  the  multiple 
scattering  nature  of  coda  and  emphasizes  the  observed  energy  flux.  Their  model  questions 
several  important  assumptions  inherent  in  the  Aki  and  Chouet  (1975)  single  backscattering 
"Bom  approximation"  model.  Central  to  these  questions  of  the  single  scattering  model  is  the 
issue  of  conversion  between  body  and  surface  wave  modes  of  propagation.  Dainty  et  al. 
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(1987)  argue  that  late  arriving  coda  (>  10  to  15  sec)  from  local  earthquakes  is  horizontally  pro¬ 
pagating  modal  or  Lg  energy.  Der  et  al.  (1984)  propose  that  a  significant  portion  of  the  Lg 
coda  is  composed  of  short-period  fundamental  Rayleigh  waves  that  concentrate  energy  near  the 
surface  and  are  multiply  scattered,  contributing  to  long  Lg  coda  on  regional  seismograms. 

In  this  report  we  conduct  simulations  of  fundamental  scattering  experiments  of  Rayleigh 
waves  incident  upon  shallow  heterogeneity.  Our  goal  is  to  determine  some  basic  quantities 
such  as  scattering  strength  and  conversion  efficiency  as  well  as  backscattering  strength  of  shal¬ 
low  heterogeneities.  These  quantities  are  needed  in  order  to  compose  models  for  coda 
observed  near  the  surface  of  the  Earth.  McLaughlin  and  Jih  (1986b)  examined  the  interaction 
of  incident  Rayleigh  waves  upon  2-D  topography  using  two-dimensional  linear  elastodynamic 
finite  difference  method  (LFD).  In  this  report,  we  explore  the  stability  of  planar  Rayleigh 
waves  incident  upon  a  2-D  shallow  heterogeneous  layer  of  various  thicknesses,  scale  lengths, 
and  rms  velocity  fluctuations  embedded  in  a  homogeneous  half  space  using  the  same  method 
of  analysis.  All  of  the  geometries  simulated  in  this  work  represent  planar  Rayleigh  waves 
incident  upon  2-D  heterogeneity. 

We  first  examine  the  problem  of  transmission  and  the  partitioning  of  energy  by  the 
scattering  process.  It  is  demonstrated  as  expected  that  larger  scale  variations  in  velocity  pro¬ 
duce  larger  attenuation.  Reflection  (backscattering)  is  observed  to  be  an  inefficient  process  for 


all  cases  examined  except  for  folded  media.  In  this  way  we  characterize  shallow  heterogeneity 
by  a  spatial  Q  for  the  transmitted  Rayleigh  waves  and  examine  the  partitioning  of  energy  into 
scattered  body  waves.  We  use  the  convened  body  wave  field  at  depth  in  the  half-space  as  a 
measure  of  the  simulated  far-field  P  and  S  coda.  Compared  to  the  strong  attenuation  of  short- 
period  fundamental  Rayleigh  waves  due  to  moderately  rough  topography,  we  find  that  shallow 
heterogeneity  generates  fewer  Rayleigh-to-body  wave  conversions,  and  hence  attenuation  is  less 
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for  shallow  heterogeneity  than  for  rough  topography.  Examination  of  the  converted  body  wave 
field  in  the  half-space  leads  us  to  conclude  that  much  of  the  SV  enrgy  will  be  trapped  in  the 
crust  as  Lg.  Conversely,  it  can  be  argued  that  incident  SV  and  P  waves  will  generate  Rayleigh 
waves  that  propagate  as  coda  on  the  surface. 

We  have  chosen  LED  modeling  out  of  convenience  and  consistency.  The  method  has  the 
advantages  that  the  solution  contains  all  conversions  (P-SV,  SV-P)  and  all  orders  of  multiple 
scattering.  We  use  the  2nd  order  formulation  of  Kelley  et  al.  (1976)  with  absorbing  boundary 
conditions  described  by  Clayton  and  Engquist  (1977)  and  1st  order  free-surface  boundary  con¬ 
ditions  for  topographic  surfaces  described  by  Jih  et  al.  (1986).  The  method  permits  examina¬ 
tion  of  fairly  general  2-D  models.  Surface  waves  are  introduced  as  initial  conditions  within  the 
haif-space  grid.  There  is  now  a  considerable  body  of  work  covering  many  aspects  of  scatter¬ 
ing  using  2-D  LFD  methods  (as  reviewed  by  by  Frankel,  1987). 
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THE  FINITE-DIFFERENCE  SCHEME 

The  2-D  linear  finite-difference  scheme  popularized  by  Kelly  et  nl  (1976)  is  used  with 
the  absorbing  boundary  conditions  (Clayton  and  Engquist,  1977;  Emerman  and  Stephen,  1983) 
on  both  the  sides  and  the  bottom  of  the  grid.  Even  with  the  absorbing  boundary  conditions, 
the  simulations  are  run  with  large  enough  grids  such  that  for  the  duration  of  the  simulations 
reflections  from  the  sides  do  not  seriously  contaminate  the  results. 

The  shallow  heterogeneous  layers  used  in  these  experiments  consist  of  3  types: 

1.  Gaussian  media  created  by  convolving  a  2-D  white  noise  field  with  2-D  Gaussian  filters 
with  specified  sizes  of  the  dominant  scatterers  ( i.e .  mean  grain  size,  a),  which  turn  out  to 
be  proportional  to  the  mean  spatial  correlation  distance. 

2.  folded  layers  with  sinusoidal  shape  of  specified  wavelength,  amplitude  and  velocity 
profile. 

3.  self-similar  media  generated  by  modulating  the  wavenumber-wavenumber  spectra  of  a 

2 

white  velocity  field  with  the  2-D  Fourier  transform,  - — — •,  of  a  special  Von  Karman 

1  +  JtV 

correlation  function. 

These  heterogeneous  media  are  then  embedded  into  a  homogeneous  half-space  with  P- 
and  S-wave  velocities  of  5.0  and  2.96  km/sec  respectively,  corresponding  to  a  Poisson’s  ratio 
of  0.35  and  a  Rayleigh  wave  velocity,  y,  of  2.71  km/sec.  Similar  random  media  have  been 
widely  used  in  modeling  the  scattering  of  either  acoustic  or  elastic  wave  by  the  finite-difference 
method  (Frankel  and  Clayton,  1984,1986;  [.evander  and  Hill,  1985;  Levander,  1985;  Frankel 
and  Wennerberg,  1987;  McLauglin  et  al.  1987a,c) 
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The  formulation  generating  the  incident  fundamental  Rayleigh  wave  packet  in  the  homo¬ 
geneous  portion  of  the  medium  is  analogous  to  that  of  Boore  (1970)  or  Munasinghe  and  Far- 
nell  (1973),  which  has  been  used  extensively  by  Martel  et  al.  (1977),  Fuyuki  and  Matsumoto 
(1980),  Fuyuki  and  Nakano  (1984),  Levander  (1985),  McLaughlin  and  Jih  (1986b).  The  initial 
waveform  is  the  fundamental  mode  Rayleigh  wave  in  the  homogeneous  half-space  with  a 
Ricker  wavelet  shape.  Figures  1  through  3  show  the  displacement  eigenfunctions  of  the  funda¬ 
mental  mode  Rayleigh  wave  (Figure  1),  the  Ricker  wavelet  and  its  wavenumber  spectrum  (Fig¬ 
ure  2),  and  the  near-surface  retrograde  particle  motion  (Figure  3).  Most  of  the  energy  of  the 
incident  wave  is  concentrated  in  the  frequency  band  below  0.9  cycle/km  =  2.4  Hz  (Figure  2). 
To  reduce  the  grid  dispersion,  we  have  oversampled  the  incident  wave  packet  by  choosing  the 
dominant  wavelength  to  be  3.2  km  =  32  Ax,  where  Ax  =0.1  km  is  the  grid  spacing.  The  whole 
grid  has  dimension  390  by  260,  with  the  heterogeneous  media  embedded  in  the  center  part  near 
the  surface. 
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NUMERICAL  EXPERIMENTS 

LFD  simulations  were  first  performed  for  24  combinations  of  shallow  Gaussian  hetero¬ 
geneities  with  4  rms  velocity  fluctuations  (5%,  7%,  10%  and  20%),  2  dominant  scatterer  sizes 
(lkm  and  2km),  and  3  different  depths  of  the  heterogeneous  layer  (1km,  2km,  and  3.2km).  10 
fixed  seismograms  are  used  to  record  the  near  surface  displacements  and  strains.  Figures  6A 
through  6H  show  snapshots  of  the  displacement  fields  due  to  a  Rayleigh  wave  propagating  in 
media  with  shallow  heterogeneity  of  fixed  thickness  h  and  various  velocity  fluctuations  v  and 
mean  scatterer  sizes  a.  Successive  frames  are  separated  by  intervals  of  2  seconds.  The 
reflection  and  transmission  coefficients  of  Rayleigh  waves  incident  upon  such  media  are  depen¬ 
dent  on  all  these  factors:  v,  a,  h.  Here  we  define  the  energy  transmission  coefficient,  T,  as 
Power, raru(f) 


T{f)  = 

P(f)  = 


P  owermc(J) 
__  P°werrtf(j) 


The  reflection  coefficient,  R,  is  defined  in  an  analogous  way: 


Power  spectra  were  used  to  compute  the  coefficients  so  as  not  to  be 


PowerMC(f) 

biased  by  the  dispersion  effects.  The  strain  in  each  case  was  "measured"  one  grid  point  below 
the  flat  surface  away  from  the  heterogeneity  in  question. 

McLaughlin  and  Jih  (1986b)  conducted  experiments  to  justify  that,  in  the  case  of  a  homo¬ 
geneous  half-space,  the  grid  dispersion  only  modulates  the  spectra  of  the  propagating  Rayleigh 
wave.  They  also  showed  that  spectral  ratios  can  be  used  to  estimate  the  scattering  attenuation 
due  to  random  topography.  We  computed  the  attenuation,  T,  for  the  shallow  heterogeneities  in 

a  similar  way  as  did  McLaughlin  and  Jih  (1986b):  V(J)  =  Q~x(f)  =  -ln(T(f)){ — % — )  where  y  = 

27T/A 

2.71  km/sec  is  the  group  velocity  of  the  incident  Rayleigh  wave  and  where  A  the  total  length 
of  the  heterogeneity,  which  is  fixed  to  be  12  km  in  all  our  experiments,  'fable  1  below  gives 
the  effective  Q  at  0.78  Hz  measured  with  vertical  displacement  near  the  surface.  The  attenua¬ 
tion  factor,  1  IQ,  plotted  as  a  function  of  frequency,  is  shown  in  Figure  8  for  these  eight  cases. 
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TABLE  1.  Q  measured  at /  =  0.78Hz 


/ 

V 

a 

Qv 

.78 

5% 

i 

617 

.78 

7% 

i 

248 

.78 

10% 

i 

127 

.78 

20% 

i 

38  ±  18% 

.78 

5% 

2 

1471 

.78 

7% 

2 

394 

.78 

10% 

2 

185 

.78 

20% 

2 

58  ±  9% 

For  larger  velocity  variation  (e  g.  10%  or  20%),  the  attenuation  tends  to  be  insensitive  to 
the  thickness  of  the  heterogeneity,  and  the  Q  measurement  would  be  more  accurate.  Frankel 
and  Clayton  (1984)  and  McLaughlin  et  al.  (1987a)  reported  that  Q~x  was  proportional  to  the 
variance  of  the  random  medium  for  body  waves  in  media  where  the  scatterers  are  comparable 
to  the  seismic  wavelength.  This  is  consistent  with  result  shown  in  Table  1  in  which  Qv  varies 
in  a  manner  essentially  proportional  to  v2. 

A  standard  (least-squares)  regression  calculation  is  performed  assuming  that  the  attenua¬ 
tion  factor  is  proportional  to  a  power  of  the  equivalent  scattering  energy-flux,  (see  Appen¬ 
dix): 

r‘ek=r"5’  (,) 

The  equivalent  scattering  energy-flux  of  the  medium,  £,  is  defined  as 

4  =  <  J  (/«2  +fJ)  dz  >  (2) 

where  the  equivalent  scattering  forces,  fu  and  fw,  are  derived  from  a  perturbation  analysis  to  the 
equations  of  motion  due  to  the  perturbations  in  the  elastic  moduli,  5(^+2|i)  and  5p, 

/„  =  8(X  +  2 p)  [Ua  +  W„Y+  8p  ((/„  -  WJ  (3a) 
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+  2|i)  Wa  +  WJ  +  6p  IW„  -  f/„]  (3b) 

and  where  (/„,  Un,Uzt,  etc,  are  the  2nd  order  spatial  derivatives  of  the  displacement  eigen¬ 
function  of  the  fundamental  Rayleigh  wave  in  a  homogeneous  half  space.  A  more  detailed 
derivation  of  the  above  formulas  is  contained  in  the  Appendix.  Under  the  assumption  that  the 
medium  has  constant  density,  p,  across  the  grid,  5(X  +  2|i)  and  5|i  reduce  to  the  velocity 
fluctuations  a2  -  a2  and  p2  -  p02  respectively. 

For  fixed  ka,  the  attenuation  factor,  1  IQ,  increases  with  the  energy-flux  of  the  medium 
(Figure  9).  But  the  rate  of  increase  depends  on  ka.  We  see  that  a  first  order  perturbation 
analysis  of  the  body-forces  should  be  sufficient  to  describe  the  attenuation  of  the  Rayleigh 
waves  for  heterogeneity  up  to  20%  rms  velocity  variation.  The  eigenfunctions  U(z)  and  W(z) 
do  not  appear  to  be  perturbed  by  the  random  heterogeneity  to  the  extent  that  the  scaling  of 
equation  (1)  breaks  down.  The  attenuation  is  therefore  dependent  upon  a  kernel  that  depends 
only  on  the  distribution  of  the  scatters  with  depth  and  the  unperturbed  eigenfunctions  as  a 
function  of  depth  (equations  2  and  3a, b). 

A  horizontal  array  of  64  sensors  are  equally  spaced  (A*  =  0.5  km)  near  the  bottom  of  the 
grid  to  record  the  converted  P  (dilatation)  and  SV  (rotation)  strain  fields.  Figures  10A  through 
10F  show  the  seismic  sections  recording  the  converted  P  wave  (dilatational  strain,  upper)  and  S 
wave  (rotational  strain,  lower)  at  every  other  sensor  of  the  horizontal  arrays.  The  coda  lasts  as 
long  as  it  takes  the  Rayleigh  wave  to  propagate  across  the  medium. 

A  frequency-wavenumber  analysis  technique  is  used  to  produce  F-K  plots  of  the  scattered 
dilatational  and  rotational  body  wavefield  as  observed  near  the  bottom  of  the  grid  (Figures  1 1 A 
through  1 1 F).  Note  that  at  3  and  4  Hz,  the  body  waves  are  mostly  forward-scattered  waves. 
The  spectra  of  the  body  waves  reflect  the  spatial  attenuation  spectra,  (?-1(f).  of  the  Rayleigh 
wave. 
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Several  observations  can  be  made  immediately.  From  0.5  to  2  Hz  the  power  level 
increases  as  h  or  v  increases.  That  those  media  with  smaller  grain  sizes  a  would  cause  more 
scattering  loss  is  also  consistent  with  the  Q  measurements,  i.e.,  Gaussian  media  with  smaller 
grain  sizes  produce  lower  Q’s.  Also,  the  converted  SV  waves  contain  substantially  more  high 
frequency  components  (2  Hz)  than  the  P  waves. 

For  Gaussian  media  with  the  same  v  and  a,  the  power  ratios  of  converted  P  to  S  waves 
appears  to  be  a  monotone  function  of  the  thickness  of  the  heterogeneity,  h.  At  1  Hz,  P/SV 
ratio  decreases  as  h  increases,  while  at  2  Hz,  P/SV  ratio  increases  with  h. 

For  comparison,  several  LFD  experiments  were  done  with  folded  layers  as  well  as  ran¬ 
dom  media  with  different  spatial  autocorrelation  structures.  Each  folded  homogeneous 
sinusoidal  layer  has  wavelength  2km,  (peak-to-peak)  amplitude  2.5km,  and  thickness  0.5km. 
The  velocity  profile  has  a  normal  distribution  with  mean  of  5  km/sec  and  standard  deviation  of 
10%.  An  interesting  phenomenon  is  obvious  from  the  snapshots  (Figure  6J),  namely  that  the 
converted  body  wavefields  are  relatively  clean  as  compared  to  all  other  cases.  Multiplely  scat¬ 
tered  Rayleigh  waves,  both  reflected  and  transmitted  ones,  are  generated  while  the  incident 
wave  propagates  across  the  layered  structure.  There  was  no  head  wave  generated  (Figure  7J), 
and  the  attenuation  measured  is  insignificantly  small.  This  is  the  only  model  which  obviously 
generates  more  back  scattered  body  waves  at  1Hz  (Figure  10J,  1 1J). 

While  for  the  Gaussian  or  exponential  media  cases,  the  correlation  distance  a  is  propor¬ 
tional  to  the  dominant  size  of  the  scatterers,  von  Karman  media  contain  scatterers  with  size 
smaller  than  a,  and  they  are  referred  to  in  Frankel  and  Clayton  (1986)  as  self-similar  media. 
A  rudimentary  comparison  indicates  that  the  thickness  of  the  self-similar  media  significantly 
affects  the  attenuation.  It  is  not  surprising  that  self-similar  media  attenuate  the  Rayleigh  wave 
more  than  Gaussian  media  of  the  same  velocity  fluctuation,  v,  and  depth. 
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DISCUSSION  AND  CONCLUSIONS 

In  order  to  use  the  results  that  we  have  generated  from  the  LFD  simulatioas,  we  need  to 
develop  an  approximate  theory  for  the  prediction  of  coda  by  the  scattering  of  Rayleigh  waves. 
In  order  to  do  this  we  will  use  the  dilatation  and  rotational  strain  sensors  located  near  the  bot¬ 
tom  of  the  grid  and  the  attenuation  estimates  for  the  Rayleigh  wave  to  give  an  estimate  of  the 
energy  radiated  by  the  scatterers  per  unit  distance  for  an  incident  Rayleigh  wave  of  unit  verti¬ 
cal  amplitude  at  the  surface.  We  will  use  a  perturbation  approach  to  replace  the  scatterers  by  a 
surface  layer  of  equivalent  sources.  The  concepts  loosely  follow  the  formalism  of  Greenfield 
(1971). 

For  our  purposes  it  is  convenient  to  think  of  the  simulation  as  though  we  had  a  continu¬ 
ous  wave  (CW)  source  radiating  Rayleigh  waves  from  the  left  side  of  the  grid  and  that  we 
observe  the  power  level  of  coda  in  the  far-field  along  the  array  of  dilatational  and  rotational 
"sensors"  near  the  bottom  of  the  grid.  The  total  energy  of  the  incident  Rayleigh  wave  packet 
in  our  simulation  is  known,  and  the  simulation  lasts  long  enough  for  the  wave  packet  to  com¬ 
pletely  clear  the  area  of  scatterers.  We  estimate  the  P  coda  power  levels  in  the  far-field  with 
respect  to  the  incident  power  by  computing  power  spectra  and  correcting  the  power  levels  for 
the  appropriate  time  windows.  The  ratio  of  the  observed  P-wave  or  S-wave  coda  power  spec¬ 
tra  to  the  input  power  spectrum  of  the  Rayleigh  wave  of  known  amplitude  then  gives  us  a 
measure  of  the  scattering  of  the  heterogeneity  as  a  function  of  frequency. 

The  dilatational  "sensors"  are  omnidirectional  detectors  of  traveling  P  waves  such  that  the 
dilatation,  0(co),  is  proportional  to  the  P-wave  amplitude,  u( co). 

02  =  (kPu)2  =  (-^-)2 

For  a  given  frequency,  the  average  P-coda  power  density,  EPc  (energy  per  unit  volume). 


Scientific  Report  #1 


12 


November  1987 


r 


AFGL-87-0322 


R-wave  Scattering  by  Heterogeneity 


TGAL-87-02 


observed  at  a  dilatational  sensor  near  the  bottom  of  the  grid  is  related  to  the  velocity  signal 
power  and  hence  ultimately  to  the  observed  dilatational  signal  power, 


EPc  -  pv2  =  pu>2u2  =  pw2*  202  =  pa202 

In  our  simulation,  the  Rayleigh  waves  were  normalized  such  that  the  Rayleigh  wave  vert¬ 
ical  amplitude  at  the  surface  at  each  frequency  was  unit  amplitude  times  the  spectral  amplitude 
of  the  Ricker  wavelet.  Therefore,  the  power  spectral  ratios  shown  in  Figures  12A  through  12E 
relate  the  average  dilatational  power  in  the  far-field  (distance  -  26*  1 05  cm)  as  a  function  of 
frequency  for  an  incident  planar  Rayleigh  wave  with  amplitude  of  105  cm,  incident  upon 
12*  105  cm  length  of  shallow  heterogeneity.  We  must  also  adjust  for  the  fact  that  the 
Rayleigh-wave  packet  is  4  seconds  long  whereas  the  coda  lasts  approximately  12  seconds  as 
the  Rayleigh-wave  packet  passes  through  the  12  km  of  heterogeneity.  The  correction  for  the 
attenuation  of  the  planar  Rayleigh  wave  in  12  km  is  of  second  order  in  our  simulations 

{e  Qv  >  0.6  for  Q  >  30,  f  ~  1  Hz,  A  =  12  km).  We  correct  the  coda  power  density  to  a  dis¬ 
tance  of  3.2  km  from  the  surface  of  the  half-space  as  if  it  were  radiated  by  a  distribution  of 

line  scatterers  along  the  free-surface.  This  correction  is  simply  a  factor  of  —  -  26/3.2  ~  8. 

rQ 

Therefore,  the  P-coda  energy  levels  radiating  from  the  bottom  of  the  layer  of  heterogeneity  at  a 
depth  of  3.2  km  in  Figures  12A-E  are  in  units  of  3.4*  10”5  erg/sec/cm2/(cm  of 
heterogeneity)/(cm2  incident  Rayleigh  wave).  Similarly,  the  S-wave  coda  power  ratios  are  in 
units  of  1 .2*  1 0-5  erglseclcrr?l(cm  of  heterogeneity )/(cm2  incident  Rayleigh  wave).  We  will  give 
these  two  quantities  the  designation,  «),  and  ^$(co).  In  our  simulations,  the  S-wave  coda 
energy  density  level  is  generally  about  equal  to  the  P-wave  coda  energy  density  level.  Since 
body  vector  forces  generate  S-wave/P-wave  amplitude  ratios  proportional  to  (a/(3)2  and  body 
stresses  generate  S-wave/P-wave  amplitude  ratios  proportional  to  (a/())\  we  see  that  the 
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equivalent  sources  required  to  describe  the  scatterers  can  not  be  simply  described  as  either 
body  forces  or  body  stresses.  The  answer  may  lie  in  the  multiple  scattering  and  the  free  surface 
conversion  of  P  to  S  and  S  to  P.  This  implies  that  we  will  require  a  separate  set  of  equivalent 
sources  for  the  P-coda  and  S-coda  until  a  self-consistent  model  can  be  constructed. 

For  example,  we  use  the  results  from  the  simulation  with  20%  rms  velocity  variation, 
Gaussian  scale  length  of  1  km,  in  a  surface  layer  3.2  km  thick.  The  P-coda  energy  level  at  1 
Hz  is  2.6*  10-6  erg!  sec!  err? I  cm!  cm 2.  This  is  equivalent  to  an  explosive  point-source  with  rms 
moment-rate  of  0.2*  1020  dyne-cm/sec  at  a  distance  of  3.2  km. 

In  order  to  predict  the  Rayleigh  to  P-coda  energy  density  we  must  predict  the  Rayleigh- 
wave  amplitude  at  a  given  distance  from  the  epicenter.  The  Rayleigh-wave  vertical  displace¬ 
ment  amplitude  from  an  explosive  point  source  at  depth,  h,  observed  at  a  distance,  A,  can  be 
written  as  (Hudson  and  Douglas,  1975) 

Uz((0)  =  A(a,P)  B(a,p,p,C0)  A- 1/2  exp(-rnhr j/y)  exp{  -  -2^- )  m(co) 


B(a,p,p,co)  = 


27tp  8pa2p3 


A(a,P)= 


2  -  f/p2  -  7171'  -  (ri'/rj  +  r|/rT)/2 


Where  a,  p,  and  y  are  the  P-,  S-,  and  Rayleigh-wave  velocities  and  where  i)  =  (l~y2/a2)1/2 
and  r|'  =  (l-y^p2)1'2.  QR  is  the  Rayleigh  wave  quality  factor.  For  example,  a  1020  dyne-cm 
explosive  source  buried  0.5  km  will  produce  3*10~4  cm  of  1  Hz  motion  at  10  km  distance  for 
our  model  half-space.  So,  our  previous  example  becomes  an  equivalent  explosive  source  den¬ 
sity  of  6*10~15  dyne-cm/sec/cm  or  6*1020  dyne-cm/sec/km. 
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We  can  now  form  a  simple  approximate  model  for  the  production  of  teleseismic  P  coda 
from  an  explosion.  As  Greenfield  (1971)  pointed  out,  the  P-coda/P  power  ratio  can  be 
predicted  as  a  function  of  time  if  the  conversion  efficiency  (scattering  cross-section)  is  known 
for  short-period  Rayleigh  to  downward  scattered  P-wave  energy.  For  teleseismic  paths,  the 
geometrical  spreading  for  P-coda  and  the  direct  P  approximately  cancel  in  the  P-coda/P  ratio. 
The  important  factors  are  the  relative  excitation  of  P-wave  amplitude  and  Rayleigh-wave 
amplitude  by  the  coda  sources  and  the  explosion  source,  respectively,  integrated  over  the 
scattering  volume  per  unit  time. 

For  the  case  of  the  Rayleigh  wave  spreading  out  from  a  point  source,  the  volume  is 
replaced  by  the  surface  area  at  time  t  =  ~  contributing  equivalent  scattering  sources;  dA  = 

2itAdA  =  lirftdt.  We  assume  that  the  energy  is  additive  from  the  scattering  surface  and  that  to 
first  order  the  scattering  is  isotropic.  We  can  modify  this  assumption  later  by  the  introduction 
of  a  radiation  pattern  for  the  equivalent  scatterers  and  the  computation  of  numerical  integrals 
over  the  area  swept  out  in  each  increment  in  time  by  the  Rayleigh  wave.  However,  to  first 
order  we  can  write  that  the  P-coda  power  received  in  a  time  interval  dt  and  frequency  interval 
dot  normalized  to  the  direct  P-wave  is  given  by 

<u2Pc>  ,  ,  -(^L>  ,  , 

— - — dtdw  =  A2(a,P)C2(a,p,p,(u)  ^co)  e  Y  (yt)  e  y  Ixy2  t  duka 

Up 

We  see  that  if  the  scatterers  are  evenly  distributed  then  the  time  decay  of  the  P-coda  gen¬ 
erated  in  this  manner  is  controlled  by  the  QR  of  the  short-period  Rayleigh,  while  the  frequency 
dependence  is  determined  by  the  depth  of  burial,  h ,  and  the  frequency  dependence  of  the 
scattering  ^(to).  This  suggests  that  if  all  P-coda  were  generated  in  this  manner,  we  could 
determine  the  Q  from  the  coda  time  decay.  Since  there  is  coda  generated  by  other  mechan¬ 
isms,  this  is  probably  not  a  very  helpful  procedure.  Also,  since  the  direct  P  spectra  may  be 
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contaminated  by  P+pP  interference  effects,  we  should  more  properly  consider  the  normalization 
factor  as  a  more  idealized  explosion  source  spectrum.  The  actual  level  of  generated  P-coda  is 
a  complicated  function  of  h,  and  the  functions  A  and  B.  For  layered  media  we  can  expect 
that  the  dependence  will  become  even  more  complicated. 

For  the  20%  rms  Gaussian  model  with  1  km  scale  length  heterogeneity  in  the  upper  3.2 
km,  we  can  predict  the  P-coda/P  signal  energy  as  a  function  of  time.  Greenfield  (1971)  com¬ 
puted  the  average  coda  signal  power  in  the  10  to  40  second  window  from  a  bandpass  filter 
normalized  to  the  total  energy  in  the  first  10  seconds  of  the  P  wave,  Wq,  given  by 

40 

j  S\t)dt 

Wn  =  _LJ° _ 

0  30  10 

His  Vk0  estimates  for  1  Hz  from  two  Novaya  Zemlya  explosions  of  21  Oct  1967  and  27  Oct 
1966  on  LASA  beams  were  about  0.0085  and  0.015  respectively.  Since  the  beams  reject  local 
coda  this  is  probably  as  reasonable  a  measure  of  the  near  source  generated  coda  as  can  be 
obtained.  Our  20%  rms  Gaussian  model  predicts  =  0.011,  which  is  in  the  right  neighbor¬ 
hood  given  the  crudeness  of  our  model.  Velocity  heterogeneity  of  20%  rms  in  the  upper  3.2 
km  probably  is  an  upper  limit  if  we  consider  typical  velocity  distributions  of  sedimentary 
rocks.  Since  there  are  other  mechanisms  of  near  source  scattering,  this  may  not  be  unreason¬ 
able  result.  We  have  found  that  rough  topography  produces  greater  Rayleigh-wave  scattering 
than  shallow  heterogeneity,  and  therefore  perhaps  topography  is  responsible  for  more  scattering 
at  Novaya  Zemlya  than  is  shallow  heterogeneity. 

For  comparison,  the  following  table  compares  predicted  Wo  values  for  several  of  our  shal¬ 
low  heterogeneity  models  at  1  Hz.  "SSI"  and  "SS2"  are  the  self-similar  models  with  1  and  2 
km  scale  lengths  respectively.  The  "GA1"  and  "GA2"  models  are  the  Gaussian  autocorrelation 
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models  with  scale  lengths  of  1  and  2  km  respectively.  The  rms  velocity  variation,  v,  is  indi¬ 
cated  for  each  model  as  well  as  the  predicted  value  of  C,p  and  the  Rayleigh  wave  Q.  We  see 
that  there  are  several  models  that  predict  values  of  W0  in  a  range  around  0.01.  Since  Q(f)  and 
C(/)  are  roughly  flat  between  0.5  and  1.5  Hz  for  most  of  the  models,  we  would  expect  that  the 
frequency  dependence  of  the  Rayleigh  to  P-coda  ratios  would  reflect  the  terms  for  excitation  of 
the  Rayleigh  wave,  to3/2e"*“to',1),  and  the  attenuation  of  the  short-period  Rayleigh  wave, 
e  Greenfield’s  LASA  beam  data  for  Novaya  Zemlya  in  the  0.5  to  2.0  Hz  range  suggest 

a  Qr  in  the  range  of  20  to  80  if  we  assume  that  the  dependence  of  W0(/)  reflects  the  depen¬ 
dence  of  QR(f).  Either  the  self-similar  models  or  the  Gaussian  models  with  rms  velocity  varia¬ 
tions  in  the  7  to  15%  range  and  scale  lengths  in  the  1  to  2  km  range  could  be  found  to  match 
his  results.  Until  such  time  as  more  detailed  analysis  is  made  of  the  P-coda/P  as  a  function  of 
frequency  and  time,  there  are  far  too  few  degrees  of  freedom  to  constrain  the  selection  of 
scattering  models.  Of  particular  importance  would  be  array  analysis  of  P-coda  using  band-pass 
filters  applied  to  multichannel  deconvolution  results  such  as  those  of  Der  et  al.  (1987).  Also, 
analysis  of  the  type  performed  by  Lay  (1987b)  to  image  scatters  using  broadband  data  coupled 
with  forward  models  would  help  to  delineate  the  actual  nature  of  the  scatterers. 
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1 

Hz  W0  PREDICTIONS 

MODEL 

V 

^0 

C,. 

Qr 

SSI 

20% 

1.25544e-03 

2.14525e-06 

19.95 

SSI 

10% 

3.66908e-02 

1.86844e-06 

70.79 

SSI 

7% 

5.43877e-02 

1.554 10e-06 

104.71 

SSI 

5% 

6.020 12e-02 

1.35356e-06 

128.82 

SS2 

20% 

1.02967e-03 

2. 14525e-06 

19.05 

SS2 

10% 

4.04064e-02 

1.70404e-06 

79.43 

SS2 

7% 

6.36782e-02 

1.62734e-06 

114.82 

SS2 

5% 

6.46347e-02 

1.35356e-06 

138.04 

GA1 

20% 

1.09775e-02 

2.579 16e-06 

35.48 

GA1 

10% 

3.75454e-02 

1.91 196e-06 

70.79 

GA1 

7% 

4.73805e-02 

1.62734e-06 

91.20 

GA1 

5% 

5.10469e-02 

1.41736e-06 

107.15 

GA2 

20% 

9.08656e-03 

1 . 1789 1 e-06 

44.67 

GA2 

10% 

2.325 15e-02 

9.80571e-07 

79.43 

GA2 

7% 

2.61012e-02 

8.15603e-07 

97.72 

GA2 

5% 

2.60344e-02 

6.47857e-07 

117.49 

The  analysis  of  the  P-coda/P  ratio  allows  us  to  cancel  the  propagation  term  for  the  P-coda 
from  scatterer  to  receiver.  In  order  to  use  the  data  for  analysis  of  the  potential  excitation  of 
SV-Lg  by  the  conversion  of  Rayleigh  to  SV,  Green’s  functions  will  be  required  for  the  SV-Lg 
wave  packet.  This  is  a  somewhat  more  difficult  problem,  but  it  can  in  theory  be  done  much 
the  same  way  as  we  have  treated  the  P-coda  problem.  This  important  problem  is  left  for  a 
later  analysis. 

We  see  that  moderate  heterogeneity  in  a  half  space  does  not  attenuate  short-period  funda¬ 
mental  Rayleigh  waves  nearly  so  much  as  rough  topography  does  (McLaughlin  and  Jih, 
1986b),  but  it  can  still  contribute  substantial  P-coda  and  moderate  attenuation  of  Rg.  For  the 
Gaussian  media  used  in  these  simulations,  the  energy  lost  due  to  body  wave  conversion  varies 
from  several  percent  to  20%  in  12  km  distance.  A  significant  result  of  the  simulations  is  that 
reflection  of  Rayleigh  waves  by  heterogeneity  at  normal  incidence  is  in  most  cases  inefficient, 
as  was  the  case  for  rough  topography.  The  only  exception  observed  was  a  folded  structure 
with  a  resonant  response  to  the  incident  Rayleigh  wave.  Therefore  we  should  not  expect  to  see 
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Rayleigh-wave  back-scattering  as  a  significant  contributor  to  the  multiple  scattering  of  funda¬ 
mental  Rayleigh  waves  that  populate  coda  for  Gaussian  or  self-similar  media.  Backscattering 
can  be  significant  for  media  with  well  defined  organized  structures  such  as  folded  sedimentary 
structures.  In  such  a  case  the  backscattered  wave  has  a  narrow  narrow  bandwidth  reflecting 
the  resonance  phenomenon. 

More  complicated  random  media  contain  many  scale  lengths  and  introduce  broad  band 
effects.  The  Rayleigh-wave  attenuation  is  a  complicated  function  of  frequency,  but  it  attains  a 
maximum  in  the  range  where  the  characteristic  wavelength  of  the  medium  matches  the 
wavelength  of  the  incident  Rayleigh  wave. 

At  low  frequency  (^_>  a )  the  coda  dilatational  and  rotational  wavenumber  spectra  indicate 
that  the  scattered  P  and  SV  waves  are  scattered  in  the  forward  direction  except  for  strongly 
folded  structures  or  Gaussian  medium  with  very  strong  velocity  fluctuations  (20%).  For  higher 
frequencies  (X  <  a  Hz),  the  scattered  body  waves  are  always  maximum  in  the  forward  direc¬ 
tion.  A  detailed  analysis  of  Rayleigh  to  P-coda  scattering  will  h  ve  to  take  this  effect  into 
account  with  an  effective  radiation  pattern  to  the  "equivalent  scatterer". 

The  results  presented  offer  a  beginning  approach  to  understanding  the  effects  upon  sur¬ 
face  waves  of  scattering  by  lateral  heterogeneity.  A  complete  exploration  of  the  problem  will 
require  variation  of  the  P  and  S  velocities,  near  surface  velocity  gradients,  and  crustal  velocity 
heterogeneity  as  a  function  of  depth. 

Much  of  the  SV  energy  scattered  by  near  surface  heterogeneities  is  concentrated  at 
apparent  velocities  within  150%  of  the  Rayleigh  phase  velocity.  A  crust  with  surface  layer  with 
P  =  2.96  km/sec  as  used  in  these  simulations  would  leak  much  of  the  energy  to  the  mantle. 
However,  if  the  near-surface  velocity  of  the  model  is  lowered,  then  the  slowness  space  occu¬ 
pied  by  the  scattered  waves  will  scale  to  the  Rayleigh  phase  velocity  and  more  energy  will  be 
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trapped  in  the  crust.  Other  ways  of  increasing  the  trapping  of  the  scattered  SV  energy  are 
increasing  the  a/p  ratio  in  the  near  surface,  introduction  of  gradients  near  the  surface  to  create 
higher  order  modes,  and  introduction  of  deeper  velocity  heterogeneity  to  scatter  P  and  SV 
energy  back  into  the  waveguide.  In  short,  all  these  mechanisms  can  act  only  to  increase  the 
Rayleigh-to-Lg  coupling.  Therfore,  we  expect  that  in  real  seismological  situations  much  of  the 
Rayleigh-wave  energy  scattered  into  SV  by  near-surface  heterogeneity  will  be  trapped  in  the 
crust  and  will  find  a  path  to  the  Lg  wavepacket. 
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FIGURE  CAPTIONS 


Figure  1  Rayleigh  wave  displacement  eigenfunctions  for  the  vertical  and  horizontal  com¬ 
ponents  of  motion. 

Figure  2A  Wavenumber  spectrum  S(/k)  =  (^^0)2exp[l-(/k//k0)2l  of  Ricker  wavelet  where  k  is 
wavenumber  and  kg  =  0.3125  cycle/km  the  dominant  (i.e.  maximum  amplitude)  wavenumber. 

Figure  2B  Ricker  wavelet  in  the  spatial  domain.  IFFT  was  applied  on  S(k)  to  speed  up  the 
wavelet  generation. 

Figure  3  The  retrograde  particle  motion  near  the  surface. 

Figure  4  Synthetic  near-surface  vertical  displacements  (upper)  and  horizontal  displacements 
(below)  for  a  Rayleigh  wave  propagating  in  a  homogeneous  half-space. 

Figure  5  Spectral  ratio  of  the  transmitted  Rayleigh  waves  over  the  incident  Rayleigh  waves  in 
a  homogeneous  half-space.  Grid  dispersion  only  modulates  the  spectra. 

Figure  6A  The  snapshots  of  the  displacement  field  due  to  Rayleigh  wave  propagating  in  a 
medium  with  shallow  heterogeneity  of  v  =  5%,  a  =  lkm,  h  =  3.2 km.  Successive  frames,  (A) 
through  (H),  are  separated  by  2  sec  intervals.  Displacements  are  proportional  to  the  darkness 
of  the  plot  and  are  normalized  to  the  maximum  in  each  frame. 

Figure  6B  Same  as  Figure  A6  except  of  v  =  5%,  a  =  2km,  h  =  1.2km.  The  diffracted 
wavefield  is  less  complicated  as  compared  to  Figure  6A  because  smaller  grain  size  causes  more 
reflection  ((C), (D)  of  Figure  6A). 

Figure  6C  Same  as  Figure  6A  except  of  v  =  7%,  a  =  1  km,  h  =  3.2 km. 

Figure  6D  Same  as  Figure  6 A  except  of  v  =  7 %,  a  =  2 km,  h  =  3.2 km. 

Figure  6E  Same  as  Figure  6 A  except  of  v  =  10%,  a  =  1  km,  h  -  3.2 km. 

Figure  6F  Same  as  Figure  6A  except  of  v  =  10%,  a  =  2 km,  h  =  3.2 km. 

Figure  6G  Same  as  Figure  6A  except  of  v  =  20%,  a  =  1  km,  h  =  1.1km. 

Figure  6H  Same  as  Figure  6A  except  of  v  =  20%,  a  =  2 km,  h  =  3.2 km. 

Figure  61  Same  as  Figure  6A  except  of  self-similar  medium  with  a  =  2 km,  h  =  3.2 km, 
v  =  10%. 

Figure  6J  Same  as  Figure  6A  except  of  folded  sinusoidal  layers  with  h  =  3.2 km,  v  =  10%, 
X  =  2 km,  peak-to-peak  amplitude  2.5km. 

Figure  7A  Synthetic  near-surface  vertical  displacements  (upper)  and  horizontal  displacements 
(below)  for  a  Rayleigh  wave  propagating  in  a  medium  with  shallow  heterogeneity  of 
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V  =  5%,  a  =  1  km,  h  =  3.2 km. 

Figure  7B  Same  as  Figure  7A  except  of  v  =  5%,  a  =  2 km,  h  =  3.2km. 

Figure  7C  Same  as  Figure  7A  except  of  v  =  7%,  a  =  1  Icm,  h  =  3.2 km. 

Figure  7D  Same  as  Figure  7A  except  of  v  =  7%,  a  =  2 km,  h  =  3.2km. 

Figure  7E  Same  as  Figure  7A  except  of  v  =  10%,  a  =  1  km,  h  =  3.2 km. 

Figure  7F  Same  as  Figure  7A  except  of  v  =  10%,  a  =  2km,  h  =  3.2km. 

Figure  7G  Same  as  Figure  7A  except  of  v  =  20%,  a  =  1  km,  h  =  3.2 km. 

Figure  7H  Same  as  Figure  7A  except  of  v  =  20%,  a  =  2£m,  A  =  3.2km. 

Figure  71  Same  as  Figure  7A  except  of  self-similar  medium  with  a  =  2&m,  h  =  3.2km, 
v  =  10%. 

Figure  7J  Same  as  Figure  7A  except  of  folded  sinusoidal  layers  with  h  =  3.2 km,  v  =  10%, 
A.  =  2 km,  peak-to-peak  amplitude  2.5km. 

Figure  8A  Attenuation  factor  observed  as  a  function  of  frequency  for  4  Gaussian  media 
(a  =  \km,  h  =  3.2 km)  used  in  the  LFD  simulations,  "x",  triangles,  and  circles  correspond 
to  v  =  20%,  10%, 7%,  and  5%  respectively. 

Figure  8B  Same  as  Figure  8A  except  of  a  -  2 km,  h  -  3.2km. 

Figure  8C  Same  as  Figure  8A  except  of  self-similar  media  with  a  =  1  km,  h  =  3.2 km. 

Figure  8D  Same  as  Figure  8A  except  of  self-similar  media  with  a  =  2 km,  h  =  3.2km. 

Figure  8E  Same  as  Figure  8A  except  of  folded  sinusoidal  layers  with  h  =  3.2 km,  X  =  2 km, 
peak-to-peak  amplitude  2.5km. 

Figure  9A  The  attenuation  factor  1/Q  at  0.78  Hz  versus  energy-flux  Observations  are 
shown  in  circles,  and  fit  to  the  curve  T  =  T^1  2842.  All  12  random  media  (with  v  -  5%,  7%, 
10%  and  20%)  have  dominant  grain  size  1  km. 

Figure  9B  Same  as  9A  except  for  grain  size  2  km.  The  smoothed  curve  is  T  =  T^1  2991. 
Figure  9C  Same  as  9 A  except  for  attenuation  measured  at  1.56  Hz.  The  smoothed  curve  is 

r  =  roe3434- 

Figure  9D  Same  as  9A  except  for  attenuation  measured  at  1.56  Hz  and  all  media  have  grain 
size  2  km.  The  smoothed  curve  is  T  =  ro^°  19944. 

Figure  10A  Seismic  sections  recording  the  converted  P  wave  (dilatational  strain,  upper)  and  S 
wave  (rotational  strain,  lower)  at  a  line  of  32  sensors  near  the  bottom  of  the  grid  spaced  1  km 
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apart  for  the  case  of  V  =  5%,  a  =  Ikm,  h  -  1.2km.  See  Figure  6A  for  snapshots. 

Figure  10B  Same  as  Figure  10A  except  of  v  =  5%,  a  =  2km,  h  =  1.2km. 

Figure  IOC  Same  as  Figure  10A  except  of  v  =  7%,  a  =  \km,  h  =  1.2km. 

Figure  10D  Same  as  Figure  10A  except  of  v  =  7%,  a  =  2 km,  h  =  1.2km. 

Figure  10E  Same  as  Figure  10A  except  of  v  =  10%,  a  =  1  km,  h  —  1.2km. 

Figure  10F  Same  as  Figure  10A  except  of  V  =  10%,  a  =  2 km,  h  =  1.2km. 

Figure  10G  Same  as  Figure  10A  except  of  v  =  20%,  a  =  1  km,  h  -  1.2km. 

Figure  10H  Same  as  Figure  10A  except  of  V  =  20%,  a  =  2km,  h  =  1.2km. 

Figure  101  Same  as  Figure  10A  except  of  self-similar  medium  with  a  =  2 km,  h  -  1.2km, 
v  =  10%. 

Figure  10J  Same  as  Figure  10A  except  of  folded  sinusoidal  layers  with  h  =  1.2km,  v  =  10%, 
X  =  2km,  peak-to-peak  amplitude  2.5km. 

Figure  11A  Log-power  horizontal  wavenumber  spectra  for  the  dilatational  (upper)  and  rota¬ 
tional  (lower)  wavefields  recorded  across  the  bottom  linear  array  for  the  case  of  v  =  5%, 
a  -  \km,  h  -  1.2km.  Power  spectra  at  1,  2,  3,  4Hz  are  normalized  to  the  maximum.  Normally 
incident  waves  have  a  horizontal  wavenumber  of  zero.  Wavenumbers  are  indicated  for  hor¬ 
izontally  propagating  P  waves.  Forward -scattered  waves  have  negative  wavenumbers,  back- 
scattered  waves  have  positive  wavenumbers.  The  dilatational  and  rotational  strain  energy  are 
largely  confined  to  P  wave  and  S  wave  slowness  across  the  array  respectively. 

Figure  11B  Same  as  Figure  1 1 A  except  of  v  =  5%,  a  =  2 km,  h  =  1.2km. 

Figure  11C  Same  as  Figure  11 A  except  of  v  =  7%,  a  =  1  km,  h  =  1.2km. 

Figure  11D  Same  as  Figure  11 A  except  of  v  =  7%,  a  =  2 km,  h  =  1.2km. 

Figure  11E  Same  as  Figure  1 1 A  except  of  v  =  10%,  a  =  \km,  h  =  1.2km. 

Figure  11F  Same  as  Figure  11A  except  of  V  =  10%,  a  =  2 km,  h  =  1.2km. 

Figure  11G  Same  as  Figure  1 1A  except  of  v  =  20%,  a  —  \km,  h  =  1.2km. 

Figure  11H  Same  as  Figure  1 1 A  except  of  v  =  20%,  a  =  2 km,  h  =  1.2km. 

Figure  111  Same  as  Figure  1 1 A  except  of  self-similar  medium  with  a  =  2km,  h  =  1.2km, 
v  =  10%. 

Figure  11J  Same  as  Figure  1 1 A  except  of  folded  sinusoidal  layers  with 
h  =1.2km,  v  =10%,  X  =2 km,  peak-to-peak  amplitude  2.5km. 
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Figure  12A  Gaussian  autocorrelation  models,  a-  1  km,  h  -  3.2  km.  Power  spectral  ratios  of 
the  dilatational  and  rotational  strain  signals  observed  24  km  deep  within  the  grid  to  the  vertical 
displacement  of  the  incident  Rayleigh  wave.  Units  are  3.4*  10-5  and  1.2*10-5  erglseclcn?l(cm 
of  heterogeneity)/(cm2  incident  Rayleigh  wave)  for  the  P-wave  and  S-wave  coda  power  density 
at  a  depth  of  3.2  km  in  the  grid.  P-wave  coda  is  lower  set  of  values,  S-wave  coda  is  upper  set 
of  values,  v  varies  from  20%  (top),  10%,  7%,  to  5%  (bottom). 

Figure  12B  Gaussian  autocorrelation  models  with  a=  2  km.  See  Figure  12A  caption. 

Figure  12C  Self-similar  autocorrelation  models  with  a=  1  km.  See  Figure  12A  caption. 

Figure  12D  Self-similar  autocorrelation  models  with  a=  2  km.  See  Figure  12A  caption. 

Figure  12E  Folded  models.  See  Figure  12A  caption. 
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APPENDIX 

We  start  with  the  2-D  equation  of  motion  for  the  vertical  component  of  motion, 

pco2(/  =  (X  +  2p)(f/«  +  W„)  +  ft(£/„  -  UX2).  A 1 

U  and  W  are  the  horizontal  and  vertical  displacement  components  of  motion  and  Ua,  Ua,  and 
UIZ  are  the  second  order  partial  derivatives  of  the  motion  with  respect  to  range,  x,  and  depth,  z. 
The  derivation  proceeds  in  a  similar  manner  for  the  horizontal  component  of  motion,  which  is 
left  to  the  reader.  The  kinetic  energy  density  in  the  vertical  component  field  is  given  by 

E  =  -yp (o2U2.  It  follows  from  the  Virial  theorem  that  for  linear  elastic  systems  the  average 

potential  energy  is  equal  to  the  average  kinetic  energy  in  the  system,  so  the  total  energy  is 
pco  U  .  We  will  equate  the  total  energy  to  the  sum  of  the  incoherent  scattered  and  coherent 
incident  field,  and  we  will  assume  that  energy  flows  from  the  coherent  to  the  incoherent  field. 

We  let  the  vertical  displacement  field  be  given  by 

U  =  U0  +  Ux>  A2 


where  U0  is  referred  to  as  the  "background"  field  and  Ux  is  the  "scattered"  field.  The  elastic 
modulus  fields  are  given  by 

(A  +  2ft)  =  (A  +  2ft)  +  5(A  +  2(i)  A3 

and 


ft  =  ft  +  5p. 

Then  to  first  order  we  have  for  the  background  field  that 

pco2t/0  =  (A.  +  2ft)((7Ql_t  +  W0a)  +  P((/0„  - 


A4 


A5 


and 
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pa>2U{  =  S(X  +  2\i)(U0xx+W0xz)  +  6p(f/0„  -  U0zz)  A6 

+  (X  +  2n)(t/,«  +  W,„)  +  ?((/,„  -  UXxz). 

The  first  term  in  (A6)  is  a  source  term  that  represents  the  "equivalent"  forces  that  give  rise  to 
the  scattered  field.  The  second  term  is  easily  recognized  as  the  right  side  of  the  propagation 
equation  from  (Al).  If  to  first  order  the  scattered  field  and  the  incident  field  are  not  correlated, 

<UxUo>  =  0, 

then  we  can  express  the  total  energy  as  the  energy  of  the  scattered  field  plus  the  energy  of  the 
incident  (or  coherent)  field, 

£  =  p®2<{/2>  =  -jpt o2Wl  +  t/?  1  A7 

Since  the  energy  in  the  scattered  field  must  arise  from  the  coherent  field,  we  see  that  the  source 
of  energy  must  be  proportional  to  the  square  of  the  source  term  in  (A6), 

fu  =  [5(X  +  2h)0/q„  +  W0xz)}2.  A8 

Similarly,  we  have  that 

fl  =  [5(X  +  2p.)(f/oxr  +  W'oc)  +  5p(U/0„  -  (/0,,)]2.  A9 

The  assumption  is  that  energy  transfer  is  unidirectional  from  the  coherent  field  to  the  scattered 

_i  $£0 

field,  so  Q - and  5£o - 6£].  Naturally,  if  the  heterogeneity  modifies  the  mode  of  the 

E0 

incident  wave  then  replacement  of  U\  by  the  eigenfunction  of  the  half-space  Rayleigh  wave  is 
not  acceptable.  Also,  if  energy  passes  from  the  scattered  field  back  to  a  forward  propagating 
Rayleigh  wave  mode  (as  may  be  the  case  for  the  folded  layers)  then  resonance  is  established 
and  this  perturbation  method  is  inadequate. 

In  order  to  compute  the  total  equivalent  scattering  energy-flux  for  our  2-D  model,  we 
integrate  over  the  depth  and  range  for  which  the  velocity  perturbations  are  non-triviai.  We 
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divide  by  the  range  in  order  to  obtain  a  range-independent  quantity.  This  assumes  spatial  sta- 
tionarity  so  that  integration  over  range  is  the  same  as  propagation  over  a  suite  of  random  media 
with  the  same  statistics. 

5  =  +  fw)dzdx  =  <j(fu2  +f*)dz>.  A 10 

The  result  that  attenuation  is  not  directly  proportional  to  2,  for  fixed  ka  implies  that  some 
multiple  scattering  is  present  and  that  a  simple  energy  flow  from  the  slowly  decaying  Rayleigh 
wave  into  body  waves  is  not  an  entirely  adequate  conceptual  model.  Also,  we  have  assumed 
that  the  weak  velocity  heterogeneity  does  not  perturb  on  average  either  the  Rayleigh-wave 
phase  velocity  nor  the  eigenfunctions. 
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